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CONVERGENCE OF FULLY DISCRETE SCHEMES FOR DIFFUSIVE 
DISPERSIVE CONSERVATION LAWS WITH 
DISCONTINUOUS COEFFICIENT 

RAJIB DUTTA, UJJWAL KOLEY, AND DEEP RAY 


Abstract. We are concerned with fully-discrete schemes for the numerical approximation of 
diffusive-dispersive hyperbolic conservation laws with a discontinuous flux function in one-space 
dimension. More precisely, we show the convergence of approximate solutions, generated by the 
scheme corresponding to vanishing diffusive-dispersive scalar conservation laws with a discon¬ 
tinuous coefficient, to the corresponding scalar conservation law with discontinuous coefficient. 
Finally, the convergence is illustrated by several examples. In particular, it is delineated that the 
limiting solutions generated by the scheme need not coincide, depending on the relation between 
diffusion and the dispersion coefficients, with the classical Kruzkov-Oleinik entropy solutions, 
but contain nonclassical undercompressive shock waves. 


1. Introduction 


In this paper, we consider a finite difference method for the vanishing diffusive-dispersive ap¬ 
proximations of scalar conservation laws with a discontinuous flux 


ut + f ik{x), [e, xGRx (0, T), 

u^{x,0) = uq{x), a; € M, 


when e > 0 tends to zero with 0 < fi{e) >-)■ 0 as e >-)■ 0. Here TZ [e, u^] is a regularization term, 

depends upon two parameters e and /i(e) referred to as the diffusion and the dispersion coefficients, 
motivated by the equations of two-phase flow in porous media, T > 0 is hxed, : M x [0, T) i—>■ K 
is the unknown scalar map, uq the initial data, fc : K i—K is a spatially varying (discontinuous) 
coefficient, and the flux function / : i—>■ K is a sufficiently smooth scalar function (see Section 2 
for the complete list of assumptions). 

Motivated by the dynamic capillary pressure [7] , we consider in this paper the simplified model 


( 1 . 2 ) 


[e, /r(e); ul^ + /i(e )7 


with a third-order mixed derivatives term including one time derivative. Here /3 ,7 > 0 are fixed 
parameters. The equation (1.1) along with (1.2) serves as a concrete model of two phase flows in 
a heterogeneous porous medium. 

Moreover, drawing preliminary motivation from phase transition dynamics, we also consider 
the following specific form of the regularization term 


(1.3) 


[£, = £/? + /r(e)7 


with ,5 ,7 > 0 fixed. 

Furthermore, for the simplicity in the exposition, we assume that the flux function has the 
following particular form 


f{k{x),u) = k{x)f{u). 


Note that the flux k{x)f{u) has a possibly discontinuous spatial dependence through the coefficient 
k, which is allowed to have jump discontinuities. 
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The scalar conservation laws with a discontinuous flux function 


(1.4) ut + f{k{x),u)^ = 0 

is a special example of this type of problems, corresponds to the case /3 = 7 = 0. A simple physical 
model corresponding to (1.4) is the Witham model of car traffic flow on a highway (consult the 
monograph by Leveque [21]), where the spatially varying coefficient k corresponds to changing 
road conditions. Several other models such as two phase flow in a heterogeneous porous medium 
that arise in petroleum industry, the modeling of the clarifier thickener unit used in waste water 
treatment plants are also corresponding to (1-4). 

Independently of the smoothness of the initial data uq and fc, solutions to (1.4) are not necessar¬ 
ily smooth due to the presence of nonlinear flux term in the equation (1.4). Thus, weak solutions 
must be sought. 


Definition 1.1 (Weak solution). A weak solution of the initial value problem (1.4) is a bounded 
measurable function it : M x [0, T) —>■ M satisfying 


(1.5) 



((fitu -f ipxk{x)f{u)) dx dt + 


ip{x, 0)Mo(a;) dx = 0, 


for all ip e x [ 0 ,T)). 

It is well known that (weak) solutions may be discontinuous and they are not uniquely deter¬ 
mined by their initial data. Consequently, an entropy condition must be imposed to single out the 
physically correct solution. If k{x) is “smooth”, a weak solution u satisfies the entropy condition 
if for all convex functions 7 : K —>■ K 


r]{u)t + {k{x)Q{u))^ + k'{x) {r]'{u)f{u) - Qiu)) < 0, in V{R x [0,T]), 


where Q : K —)• K is defined by Q'{u) = r]'{u)f'{u). 

By standard limiting argument, this implies the Kruzkov-type entropy condition 

\u - cjj -b sign {u - c) {k{x){f{u) - /(c)))^ -b sign {u - c) f{c)k'{x) < 0, in V{R x [0, Tj), 
holds for all c S K. 

However, the notion of entropy solution described above breaks down when k{x) is discon¬ 
tinuous. In view of [14], we use the following notion of entropy solution for (one-dimensional) 
conservation laws with discontinuous flux equations with coefficients that are only spatially de¬ 
pendent. We assume that the spatially varying coefficient k{x) is piecewise with finitely many 
jumps (in k and k'), located at ■ • • ,^m- 


Definition 1.2 (Entropy solution). A weak solution u of the initial value problem (1.4) is called 
an entropy solution, if the following Kruzkov-type entropy inequality holds for all c € M. and all 
test functions 0 < "0 € I?(IR x [0, T]). 

( 1 . 6 ) 



c| ift + sign {u 


c)k{x){f{u) - f{c))iix) dxdt+ [ |iio 

Jm 


c\ ipix, 0) dx 


f fT M T 

/ / sign{u-c)k\x) f{c)ijdxdt+'^ \f{c){k+- k-)\i;{^rn,t) dt > 0. 

Jo Jo 


The last couple of decades have witnessed remarkable advances in the studies of conservation 
laws with discontinuous flux function. However, we will not be able to discuss the whole literature 
here, but only refer to the parts that are pertinent to the current paper. In case of “smooth” k{x), 
the notion of entropy solution was introduced independently by Kruzkov [18] and Vol’pert [31] 
(the latter author considered the smaller BV class). These authors also proved general existence, 
uniqueness, and stability results for the entropy solution, see also Oleinik [24] for similar results 
in the convex case fuu R 0 . 
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1.1. Diffusive Dispersive Approximation. It is well known that the conservation law (1.4) is 
derived by neglecting underlying small scale effects such as diffusion, dispersion, capillarity etc., 
and may admit physically relevant discontinuous solutions containing shock waves (non-classical 
shock) that may depend on underlying small-scale mechanisms. It has been successively recognized 
that a standard entropy inequality (due to Kruzkov, Oleinik, and others) does not suffice to single 
out such a physically relevant solution, and it is important to incorporate these small-scale effects 
in the entropy condition. In other words, additional admissibility criteria (a kinetic relation) are 
required in order to characterize these small-scale dependent non classical shock waves uniquely. 
In [9], the authors developed a framework for the existence and uniqueness of the non-classical 
shock waves that arise as limits of diffusive-dispersive approximations. 

Noting that the solutions to (1.4) can be different due to their explicit dependence on the 
underlying small scale effects, we focus on a concrete model of two phase flow in porous medium 
(for a brief derivation of this model consult [4]). The relevant small scale effect is a dynamic 
capillary pressure term, that was introduced in [7]. Compared to the standard capillary pressure 
models [1], the addition of the new term resulted in a model that contain higher-order mixed 
spatio-temporal derivatives (cf. (1.2)). 

The diffusive-dispersive model has a long tradition, starting with the analysis of linear diffusion- 
dispersion model (1.1). A pioneering study of the effect of vanishing diffusion and dispersion terms 
in scalar conservation laws, with ai-independent flux function, can be found in Schonbek [26]. The 
technique of compensated compactness was used to prove convergence toward weak solutions. 
Kondo and LeFloch [17] studied zero diffusion-dispersion limits for x-independent fluxes under an 
optimal balance between the sizes of the diffusion and dispersion parameters. LeFloch and Natalini 
[19] used the concept of measure-valued solution and established convergence results assuming that 
the diffusion dominates the dispersion. Subsequently, the approach of kinetic decomposition and 
velocity averaging [25] was introduced by Hwang and Tzavaras [11] to analyze singular limits 
including nonclassical shock waves. Furthermore, we also mention related works by Wu [32] and 
Jacobs, McKinney, and Shearer [12] which provides the first existence result of undercompressive 
shocks for the modified Korteweg-de Vries-Burgers equation. 

It is well known that the relative scaling between e and fi{e) determines the limiting behavior 
of solutions, and we can distinguish between three cases: 

• Diffusion-dominant regime ii{e) « e^: The qualitative behavior of solution is same as 
the solution of the conservation laws. 

• Dispersion-dominant regime ^(e) >> £^: In this case, high oscillations develop (as £ I 0) 
especially in regions of steep gradients of the solutions and only weak convergence is 
observed. 

• Balanced diffusion-dispersion regime: This typically corresponds to the scenario where 
/r(£) = 0{£^). Only mild oscillations are observed near shocks, and the limit solution is 
a weak solution to conservation laws. Most importantly, in this case, the solution exhibit 
non-classical behavior, as they contain undercompressive shocks. However, for the regime 
/r(£) = C’(£^), the limit solution coincides with the entropy solution determined by Kruzkov 
theory [18]. 


1.2. Numerical Schemes. It is well known that standard finite difference, finite volume and 
finite element methods have been very successful in computing solutions to hyperbolic conserva¬ 
tion laws with discontinuous coefficients, including those containing shock waves. However, we 
mention that most of these well-established methods are proven to be not good enough to capture 
nonclassical shock wave solutions numerically. This well known phenomena has been explained 
by many anthers (Hou and LeFloch [10], Hayes and LeFloch [8], and others) in terms of the 
equivalent equation associated with discrete schemes through a formal Taylor expansion. The key 
idea behind capturing nonclassical shocks is to design finite difference schemes whose equivalent 
equation matches, both, the diffusive and the dispersive terms (cf. (1-3)) in the underlying model. 
However, these schemes fail to approximate nonclassical solutions with large amplitude, especially 
strong shocks due to lack of control on higher order error terms present in equivalent equation. 
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A recent work by Ernest et al. [6] has overcome such problems by dominating higher order error 
terms in amplitude by the leading order terms of the equivalent equation. 

In another paper by Chalons and Lefloch [2] , the authors introduced a fully -discrete scheme for 
the numerical approximation of diffusive-dispersive hyperbolic conservation laws (cf. (1.1)-(1.3)) 
in one-space dimension. An important feature of their scheme is that it satisfies a cell entropy 
inequality and, as a consequence, the space integral of the entropy is a decreasing function of time. 
Moreover, they showed that the limiting solutions generated by the scheme contains nonclassical 
undercompressive shock waves. 

On the other hand, there is a sparsity of efficient numerical schemes for (1.1)-(1.2) available 
in the literature. In fact, to the best of our knowledge, this is the first systematic attempt to 
construct a provably convergent numerical scheme for (1.1)-(1.2). Having said this, there are some 
numerical experiments available in the final section of the recent paper by Coclite et al. [4] without 
rigorous analysis of the scheme. 

1.3. Scope and Outline of the Paper. In view of the above discussion, it is fair to claim 
that there are no robust and provably stable numerical schemes currently available to simulate 
the vanishing capillarity approximations of scalar conservation laws equation (1.1)-(1.2). In this 
context, we consider a fully-discrete (in both space and time) finite difference scheme for (1.1)- 
(1.2) which is provably convergent and able to capture non classical shocks quite well. Since 
diffusion-dispersion model for the conservation laws with discontinuous flux has not been studied in 
detail, we analyze a fully-discrete scheme for (1.1)-(1.3) as well. While there are several numerical 
methods which perform well in practice, perhaps better than the one presented here, (see [20] for 
a recent comparison of diferent numerical methods) we emphasize that we prove the convergence 
of the schemes proposed in this paper. Here, we mention that a detailed analysis of the scheme 
introduced by Ernest et al. [6] is beyond the scope of this paper, and will be the topic of an 
upcoming paper. 

To sum up, the schemes in the present paper have the following properties: 

(a) Approximate solutions for (1.1)-(1.2), generated by the scheme (3.2), converge to the 
unique entropy solution of (1.4) as long as fi{Ax) = o{Ax). A scheme (cf. (7.2)) has been 
formulated for (1.1)-(1.3) and the same techniques can be applied, mutatis mutandis,to 
prove convergence of approximate solutions to the unique entropy solution of (1.4). 

(b) Approximate solutions for (1.1)-(1.2) have been shown to converge to weak solutions of 
(1.4), when ii{Ax) = 0{Ax^). Moreover, we show numerically that the limiting solutions 
generated by the schemes (3.2) and (7.2) contain nonclassical undercompressive shock 
waves. 

The rest of the paper is organized as follows: In section 2, we present the mathematical frame¬ 
work used in this paper. In particular, we have used a compensated compactness result in the 
spirit of Tartar [27] but the proof is based on div-curl lemma and does not rely on the Young 
measure. Section 3 introduces the fully-discrete finite difference scheme for (1.1)-(1.2) . In sec¬ 
tion 4, we derive a priori estimates for the approximate solutions and a detailed convergence 
analysis towards weak solutions of (1.4) has been discussed in section 5. Convergence towards the 
unique entropy solution has been considered in section 6, while a brief discussion on the results for 
diffusive-dispersive approximation (1.1)-(1.3) has been addressed in section 7. Finally, numerical 
results are presented in section 8 to illustrate the performance of the designed schemes. 


2. Mathematical Framework 

In this section, we list all the assumptions on the data for the problem (1.1), and present relevant 
mathematical tools to be used in the subsequent analysis. Throughout this paper we use the letters 
C, K etc. to denote various generic constants independent of approximation parameters, which 
may change line to line, but the notation is kept unchanged so long as it does not impact the 
central idea. 

The basic assumptions on the data of the problem (1.1) are as follows: 
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A.l For the initial function wq : M i—>■ K, we assume that 

uq e n L°°(K), a < uq{x) < b, for a.e a; S M; 

A.2 For the discontinuous coefficient fc : K i—)■ K, we assume that 

k G L°°{R) n -BV[oc(K), a < k{x) < /3, for a.e a; G M; 

A.3 Regarding the flux function / : [a, /3] x [a, b] i-A K, we assume that 

u I—>■ f{k, u) G 6]), for all k G [a, /3], 

k i-G f{k, u) G C^{[a, /3]), for all u G [a, 6]; 

A.4 Furthermore, we assume that u i-G f{k,u) is genuinely nonlinear a.e. in K x [0,T], i.e., 
fuu{k{x),u) ^ 0, for a.e. u G [a,b]. 

Remark 2.1. It is worth mentioning that the Assumption A.4 is typically required in the com¬ 
pensated compactness framework. This condition also imposes a condition on the coefficient k(x). 
In fact, it implies that f{u) is genuinely nonlinear (i.e., f 0) and |fc(a;)| ^ 0, for a.e. x G K. 


Next, we recapitulate the results required from the compensated compactness method due to 
Murat and Tartar [23, 27]. For a nice overview of applications of the compensated compactness 
method to hyperbolic conservation laws, we refer to Chen [3]. Let A4(K) denote the space of 
bounded Radon measures on M and 


Co 


= € C(M) lim ijj{x) = 0 > . 

|x |—>00 


If ^ G A4(K), then 

= / Ip dp., for all -ip G Co{M.). 

JR 

Recall that p G A4(K) if and only if |(^,'0)1 < C ||0||j^oo(r) for all 0 G Co]®). We define the norm 

I1mI1a4(k) := sup| |(ai, 0)| : 0 G CoW, ||0||ioo(R) < l}. 

The space ||•||^(R)^ is a Banach space and it is isometrically isomorphic to the dual space 

of I^Co(K), II'll Loo Furthermore, we define the space of probablity measures 

Prob(K) := G A4(K) : ^is nonnegative and ||m||^(r) = l|. 

Before we state the compensated compactness theorem, we shall recall the celebrated div-curl 
lemma [23]. 

Lemma 2.1 (div-curl lemma). Let fl be a bounded open subset o/K^. Suppose 
1 . —1 2 . —2 1 . —1 I 2 . —2 

in ifffil) as Ax 0 0. Furthermore, assume that the two sequences ^ div 1 and 

I J Aa;>0 

|curl (uaj,, ^Aa;) I ® (common) compact subset of where div (ua,j,, Ua^,) = 

curl (uAxj'ci.a:) = ~ Then along a subsequence 


{uax^ ^L) • i^Ax^^Ax) ^ , in Vfn), as Ax 0 0. 

Suitably modified for our purpose, we shall use the following compensated compactness result. 
For a proof, we refer to the paper by Kenneth and Towers [13, Lemma 3.2]. 
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Theorem 2.1. Assume that A.2, A. 3 and A. 4 hold. Let C M x [0,T] be a bounded open set, 
and assume that {uax} a sequence of uniformly bounded functions such that |uAa:| < M, for all 
Ax. Set 

im{s),Qi{k,s)) = (s - cj{k,s) - f{k,c )), 

^ ^ iV 2 ik,s),q 2 {k,s)) = (j{k,s) - f{k,c), ifsik,0)f d9 

where c is an arbitrary constant. If the two sequences 

{ril{uAx)t + gi(fc(a;),UAx)x}Aa;>0 > ,UAx)t + g2(fc(a;),UAx)a;}Ax>0 

belong to a compact subset of then there exists a subsequence of {uax}ax>o converges 

a.e. to a function u € L°°(fl). 

We remark that, a feature of the compensated compactness result above is that it avoids the use 
of the Young measure by following an approach developed by Chen and Lu [3, 22] for the standard 
scalar conservation law. This is preferable as the fundamental theorem of Young measures applies 
most easily to functions that are continuous in all variables. 

The following compactness interpolation result (known as Murat’s lemma [23]) is useful in 
obtaining the Ilif^ compactness needed in Theorem 2.1. 

Lemma 2.2. Let fl be a bounded open subset o/K^. Suppose that the sequence {>CAa;}Aa;>o 
distributions is bounded in Suppose also that 

k^Ax — kli^Ax T ^2,Ax', 

where {/1i.ax}ax>o ® compact subset of (LI) and {C2 ,Ax}ax>q ® bounded subset of 

.^ioc(^^)- Then {>CAs}Aa;>o ® compact subset of 


3. A Fully-Discrete Finite Difference Scheme 


We begin by introducing some notation needed to define the fully-discrete finite difference 
scheme. Throughout this paper, we reserve Ax, At to denote small positive numbers that represent 
the spatial and temporal discretizations parameter of the numerical scheme. Given Ax > 0, we 
set Xj = jAx for j G Z, to denote the spatial mesh points. Similarly, we set = nAt for 
n = 0,l,-- - ,N, where NAt = T for some fixed time horizon T > 0. Moreover, for any function 
u = u{x,t) admitting point values, we write u" = u{xj,V^). Furthermore, let us introduce the 
spatial and spatial-temporal grid cells 

Ij = [xj_i/2,Xj+i/2), /j" = [xj-i/2,Xj+if2) X [i”,r+^). 


where Xj±if 2 = XjzLAxl2. Let D± denote the discrete forward and backward differences in space, 
i.e., 


D±Uj = ± 


Uj 1 Uj 

Ax 


The discrete Leibnitz rule is given by 


D±{ujVj) = UjD±Vj + Vj±iD±Uj 

while the summation-by-parts formula is given by 

jez jez 

Furthermore, for any function /, using the Taylor expansion on the sequence f{uj) we obtain 

/\ rp 

D±f{uj) = f'{uj)D±Uj ± 

for some between Uj±i and Uj. In other words, the discrete chain rule is accurate up to an 
error term of order Ax{D±Uj)'^. 
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Finally, let denote the discrete forward and backward difference operator in the time, i.e.. 


, <=^1 - < 
n* 7/" — _F 

^ At 


The following identity is readily verified: 

(3.1) 


n r^t n 1 r-it / n\2 / /-it n\2 

UjD^Uj =- . 


We propose the following fully-discrete (in space and time) finite difference scheme approximating 
the limiting solutions generated by the equation ( 1 . 1 )-( 1 . 2 ) 


(3.2) Dlu] + = l3AxD+D_u] + -ffi{Ax)DXD+D_i 

1 


j ez,ne No, 


(3.3) 


= 

^ Ax 


— / Uoio)d0, j e z. 


where / 3,7 > 0 are fixed parameters, and fi{Ax) i—>■ 0 as Ax i—>■ 0. More specifically, we will 
either use fi{Ax) = 0{Ax^) or ^{Ax) = o{Ax^) depending on the quest for the convergence of 
approximate solution mat; towards a weak solution or the entropy solution, respectively. 

Remark 3.1. Here we used the notation 0{Ax) to denote quantities that depend on Ax and are 
bounded above by CAx, where C is a constant independent of Ax. Likewise, we used the notation 
p{Ax) = o{Ax) to denote quantities that depend on Ax and are bounded above by C'Aa;“, where 
C is a constant independent of Ax and a > 1. 

The numerical flux corresponding to the flux function k{x)f{u) is given by 






with fc., 1 = . 
■^+2 Ax 


k{x) dx, 


where fJ'_^_l ■= /”(%i’^i-i-i) is based on a two-point monotone numerical flux, i.e., non-decreasing 
with respect to the first argument and non-increasing with respect to the second argument, and 
consistent with the actual flux, i.e., f(u,u) = f{u). Moreover, in order to maintain monotonicity 
of the scheme (3.2) without the higher order terms (corresponds to = 7 = 0), the arguments of 
the numerical flux are transposed when the coefficient k is negative. More specifically, we choose 


rn _ 

•'7+3 “ 


/(m”,m”+i), iffcj+i >0 


Summing up, the numerical flux h'^.i is given by 


— 


kj+i /(u?,u" 


^3 ' 
(uL 


t-ti 


.C(S'+i’ 


), 


if > 0 

if < 0 . 


In particular, we focus on Engquist-Osher (EO) numerical flux given by 
(3.4) f{u,v) = l:{f{u)+f{v))-l-[ \f{s)\ds. 


Remark 3.2. We have chosen to analyse the scheme (3.2)-(3.3) with EO flux because of its 
apparent simplicity. One can, however, adopt the method of proof developed in this paper and 
obtain similar results for other schemes (e.g., all monotone schemes). 

To this end, observe that the EO flux given by (3.4) is Lipschitz continuous . In fact, for f G C^, 
it has continuous partial derivatives satisfying 


(3.5) 


f-{v) = fv{v,u) < 0 < fu{v,u) = /+(m). 


using the conventional notations that a_ = min(a, 0) and a+ = max(a,0). It is also clear that 
11/'I I serves as a Lipschitz constant for EO flux. 
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For a given initial data Uq, we define the initial grid function by (3.3). Moreover, for 

the sequence | u? |. „ » ; associate the function u\x defined by 

U\x{x,t)= ^ 

jGZ,n>0 

where 1a denotes the characteristic function of the set A. Similarly, for the coefficient k approxi¬ 
mated at each cell boundary, we associate the function k/\x defined by 

kAxix) = ^3+h 

tez 

Note that k and u are discretized on grids that are staggered with respect to each other. This 
indeed results in a reduction in complexity, compared with the approach where two discretizations 
are aligned. 

Throughout this paper, we use the notation uax to denote the functions associated with the 
sequence xieNo' later use, recall that the discrete £^(IR) and norms, and 

BV semi-norm for a lattice function uax are defined respectively as 

l|uAx(',t”)||^oo(R) =SUp|u”|, ||MA:r(',t”)||fl(R) =Ax^\u]\, 

\\uAx{‘,t )I|£2(-r) = . |'a?| 1 \uAx{‘,t )Ibv ~ K-i-H ~ I 

y jGZ jez, 

For the sake of simplified notations, unless specified, we shall use the notation ||’|| to denote the 
discrete ^^(M) norm. 


4. A Priori Estimates 


This section is devoted to the derivation of a priori estimates which turns out to be useful to 
prove “strong compactness” of the approximate solution uax- To begin with, following Coclite et 
al. [4], we assume that the approximate solutions generated by the scheme are uniformly bounded, 
i.e., uax € X [0,T]). In other words, we assume that 

B.l For almost every (a;, t) G K x [0, T], a' < uax < for some fixed constants o', b' G K; 

Remark 4.1. It is worth mentioning that the above assumption on the approximate solutions is the 
manifestation of the specific structure of the flux function (depends explicitly on space variable). In 
fact, to obtain bound on the solution, one requires a priori L°° bound on the solution (cf. (4.4 ) ). 
This assumption can be toppled by replacing the “space dependent flux function” to a flux function 
which depends explicitly only on the solution. In such a scenario, one can use framework of 
the compensated compactness result [5], to reproduce all the results in this paper. 

To proceed further, we first collect all the available estimates on the approximate solutions in 
the following lemma. 


Lemma 4.1. Let uax be a sequence of approximations generated by the scheme (3.2). Moreover, 

assume that the initial data Ug lies in L^(M). Then the following estimate holds 

(4.1) 


Inf Ik" 


2 f 7 m (^2;) fAx 






DlWD.u^f + SAx\\D_u^ 


-f SAx + S-fAxp{Ax) \\DlD_u^\\ < C, 


provided At and Ax satisfies the following CFL condition 

(4.2) max|2max|||/||^,||/'||^,||fc||^|-H^, ^ (^ + jffiAx) } | ^ ™ - 5), 

with S G (0, min (1,/?)). Here A = At/Ax and the constant C > 0 is independent of Ax. 
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In particular, the estimate (4.1) guarantees following space-time estimates: 


(4.3a) 

VnGN, AxY,iu]f<C, 

j 

(4.3b) 

Ax^At'£'E^D_u]f<C, 

j n 

(4.3c) 

Ax^At'£'E^DXu]f<C, 

j n 

(4.3d) 

AtAx^ti{Ax)J2J2^DlD_u] f < C. 


j n 


Remark 4.2. In light of the CFL condition (4.2), we want to emphasize that if fi{Ax) = 
O{Ax ) (required to prove convergence towards a weak solution, cf. Theorem 5.1), then we need 
C ■= fixed. On the other hand, if p{Ax) = o{Ax'^) (required to prove 

convergence towards the entropy solution, cf. Theorem 5.1), then we need C = (Ai)i+° ) with 
a > 5, to he kept fixed. To sum up, we need a stronger CFL condition to prove convergence of 
approximate solutions towards the unique entropy solution of (1.1). To this end, we mention that 
in the subsequent analysis the CFL condition (4.2) is always assumed to hold. 

Proof. To start with, we multiply the scheme (3.2) by Axuf and subsequently sum over j G Z. 
Then, using summation-by-parts formula and the identity (3.1), we obtain 

j 3 3 


= —f5Ax''^Ax\D_Uj'^ — 'yp{Ax)Ax'y^D-ufD*^_{D-uf). 

3 3 

Note that, the identity (3.1) also implies that 

3 3 3 

Next, we move on to estimate the term T/s.x{f)- Using summation-by-parts formula, we obtain 

-Iix(/) = -i: A,t»“DA”_. = 

3 3 

= - (F(«j)-f(«"-.))] + (f(»")-u«"-i)) 

3 

= [("? -- ('"(“?) -'"("?-.))] -E {'‘m - ^-i). 


£j.Af) 

where F is the primitive of /, i.e., F' = f. A first order Taylor’s expansion together with the 
monotonicity of the numerical flux function f{u(j,u'j^i) gives us an estimate of £j^n{f)- To see 
this, notice that 

GM) = («" -- /(«”_.)), 

where it"_i lies in between uf and u^,. To proceed further, we consider the following two cases: 

J 2 J J 

Case 1: Assume that kj_i > 0, then by the definition of numerical flux /"_i = f{Uj_i,u(- ). If 
m" < u" 1 < u" 1 , then 

3 — 3—i ~ 3 — 


/"_i > /«_!,<) > 
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On the other hand, if < u"_i < u", then 

J J 2 

J 2 J 2 -J 2 J 2 J 2 

Thus, in any case we have £j,nif) < 0. 

Case 2: Assume that fc _i < 0, i.e., /"_i = If m" < u"_i < n?-!, then 

J 2 3 2 J J J J 2 ^ 

J 2 J 2 J 2 J 2 'J 2 

On the other hand, if u"_i < u"_i < u", then 

/”_i > />”_!, «7-l) > /(W”_i,w"_i) = 

J 2 J 2 J 2 J 2 “'2 

Therefore, in this case also, we have £j,n{f) < 0. Having this in mind and making use of the 
Assumption B.l, we conclude that 


(4.4) 






< BV(/c) max II < C, 


where C is a constant independent of Ax. Finally, combining all the above estimates, we obtain 


1 


\Di\\u^r- 


At 


Wlu^ + 


2 7/i(Ax) 


(4.5) 


DX\\D_u^^ 

7/i(Ax) At 


\dID_u^\\ + pAx ||T»-M”|r < C. 


Next, we multiply the scheme (3.2) by Ax and sum over j S Z to obtain, 

(4.6) 

Ax^Y.^D\u-f + Ax^Y.D.h-H^%u- 

3 3 

(4.7) = -P{Axf D_u]DX{D_u]) - 7/r(Ax)Ax2 Y ■ 

3 3 

Again, use of the identity (3.1) reveals that 

PiAxf Y D-U-DX{D.u-) = ^Y Dl{D-u]r - Y • 

3 3 3 

Next, considering the term which involves flux function, we see that 

ax^Y^-^hd>7 = Diu- 


+ “i 


= Ax 


i:( 




)/7.D7" + A»,j:t,_s(/ 




3 3 

Recall that we have chosen to work with a specific monotone flux, i.e., Engquist Osher flux. Since 
EO flux is Lipschitz continuous with Lipschitz constant ||/^||oo (c.f. (3.5)), we conclude that 


sup sup 

3 ^ 


/■ 


i+5 


< 


— M ./ M OO 


, and 


£71 £71 

Jj+l - Ji-i 


<ll/'lloo(h”-«”-l| + 


An application of Young’s inequality: For any two real numbers a and b and for every e > 0 

. 2 
ab < ea + — 


4e 


leads to the estimates 


|Dy.“l 


< eAx I 




4e 


— \\D+kf 
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and 


- f"-i) i ll/'IU II^ILi: (|0-“;‘l + |0-“?+il) \D 


* 7;” I 
+^J I 


<2e^Ax\\f\\^\\k\\^Dlu^ + 


2 , ll/'lloollfciloo 


2ei 


iD-u" 


For notational simplification, we define M := max| ||/||oo , ||/^||oo > Halloo }• Combining all above 
estimates, we arrive 


Ax + 'yAxfi{Ax) 


(4.8) 


/3Ax ..nii2 ^ Hot „,n||2 , ^Aa:,,^ ,_„2 


\DlD_u^\\^ < eAxM\\Diu^\\^ + 


4e 


\D+kr 


+ 2eiAxM^\\Dlu^\\ + 


t _ „||2 , M^Ax „ ^ _ „„2 


2ei 


D_u” 


Finally, adding (4.5) and (4.8) yields 
(4.9) 

1 


-Di ii„”f 


+ Ax{l- Ms- 2M^ei - -J \\Dlu^ 

+ ,A.,(Ax)ll-^-|^)||77i77_„"ir<C, 


where A = Al. 


We must now use the CFL condition (4.2) to conclude that for some S S (0,1) 
(4.10) 


1 /3 - ^ and 1 - Me - 2Mhi - ^ > S. 

2 2 jfi{Ax) 2 ei 2 


Note that 13 > S must hold. Furthermore, choosing e = l,ei = 0.5, the third contraint in (4.10) 
can be written as 

M + M^ + ^<l-S 


which would require M < 1. Assuming M is small enough (this can be done upto rescaling) to 

A 

2 


ensure the existence of a d S (0,1), such that 2M + | < min(l — 6,13 — 6). This in turn would 


imply 


A 


A 


M^ + M+ -<2M+-<l-S, 

I3-M^>P-2M >5 

In other words, we get the last two conditions of (4.10). Thus, the CFL condition (4.2) ensures 
all the conditions of (4.10) are satisfied. Consequently, we see that estimate (4.1) holds. 

In order to prove estimates (4.3a), (4.3b), (4.3c), and (4.3d), we multiply the inequality (4.9) 
by At and subsequently sum over all n = 0,1, • • • , A — 1 to reach 


1 II ivi|2 , f 7M(Aa;) ^ pAx^ 

2 11“ II + I 2 +2 


+ dAzAt^ \\D-i 


+ 6AxAt^^ + 6^Axfj.{Ax)At'^2 ||Zl(i_Zl_u"||^ < J^{uo), 


where 


^ 1 II ||2 , (lt^{Ax) PAx^\ 2 

■^Wo) := 2 ll^oll +( -^-+ ^— ) llH-Uol 
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1 

< - 
- 2 



This essentially finishes the proof of the lemma. 


7^(Aa;) /3Aa;^\ 

2 + 2 J 



< C. 


□ 


5. Convergence Analysis 


Having obtained all the necessary a priori bounds in the previous section, we are ready to 
prove that the approximate solutions generated by the scheme (3.2) converge strongly to a weak 
solution of (1.4), at least along a subsequence. The general strategy of the convergence proof is 
in the spirit of the one used by DiPerna [5], and it has been used in various contexts by several 
different authors. However, as we pointed out earlier, we shall use a simplified compensated 
compactness result in the spirit of [13, Lemma 3.2]. In what follows, using a priori estimates 
derived in Section 4, we first demonstrate the desired compactness of {uAx} ax>o- Then, an 
application of the compensated compactness Theorem 2.1 gives the desired strong convergence in 
Tfoc p < oo. To achieve our goal, we start with the following crucial lemma: 

Lemma 5.1 compactness). The sequence 

{?7j(A:(x),UAa:)t + %(^(a::), UAx)x}Aa;>o compact in x [O.T]), 

where pi and qi are given by (2.1), for i = 1,2. 


Proof. To begin with, let us assume that p = pi, for z = 1 or * = 2, and p he a, test function with 
compact support such that ip{x) = 0, for all |a:| > |a;j_|_i |, and for t > NAt, for some J and N. 
With this ip, we define the following functional 

i^Axip, q),T) ■= -{p{k(x), UAx)t + q{k{x), UAx)x, p) 

= / 'niHx),UAx)Tt + q{k{x),UAx)Txdxdt:= {£Ax{p,q),T) + {£Axiv,Q),T), 


where 


{^Axi'n,q),‘f) = / / {r]{H^),UAx) - p{kAx,UAx))Ttdxdt 

Jo Jr 


{q{k{x), UAx) - q{kAx, UAx)) Tx dxdt, 


and 


{SAx{p,q),'p) = / 'n{kAx,UAx)Tt + qikAx,UAx)>Pxdxdt. 

Jo Jr 

In what follows, we let Ht- denote an arbitrary but fixed bounded open subset of K x [0,T]. Let 
r G (1, 2] and set p = G [1, oo). With p G we have by Holder’s inequality 



\{£ax{v,q),t)\ < C\\k 

kAx\\LP(nT) Wo'''(nT) 

so that 




(5.1) 


{£AxiV,Q)}Ax>0 

is compact in 

Next we focus 

on 

the other term 




N-l t" + l j; 

3-¥ ^ 

{£Ax{p,q), 


"Si:/ / 

^ p(k^\vj;)pt + q{k^^, 


j n—0 * 

N-l t"+i 

= / 

j n=0'^*" 


EE 

j n=0 


p{kj_i,u'^)pt + q{kj_i,u^)px dxdt 

p{k^j^i,vJ])pt P q{kj+i,u'^)pa: dxdt 


N-l 

I 9. 
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where 

N-l .t"' 

(5.2) 

{£Axiv,q),‘p) = E E / 


j n=0 

and 

N-l (»+’ 

(5.3) 

{^Ax{V,q),<f) = E E / 


3-2 


f 3^ 2 


tpt dxdt 

N-l px 1 

-^EE r~ f(7(fcj+i,u”) - q{kj_i,vJ^)\ ipx dxdt. 

, n=o 3t" 3x^ 


The proof is essentially complete if we assume that both and {'S’ax(^> 9)}Aa:>o 

are compact in □ 

In the proof Lemma 5.1, we claimed the compactness of 9 )}ax>o '^)}Aa;>o' 

We try to justify this claim below. For the rest of this section, we assume that the conditions 
stated in Lemma 5.1 hold. We will also continue to use IIt’ introduced in the proof of the above 
lemma. 

Lemma 5.2. 

(5.4) {^L(^>9)}a2,>o is compact in H-l{nT). 

Proof. Consider the expression given by (5.3), which is split as 

i^Axiv,<]),¥’) ■= {£Alid,q),p) + {£Alid,q),v) 

where 

AT 1 

-r]{kj_i,u])^ ifit dxdt 


N-l ^ 

/ 3+i5 


{£^/xi'n,<i),p) =51 E 

j n=0 

^“1 px 1 

(4’f(.-.),^)=EE r~ fg(fcj+i,M") - q{kj_i,u^)\ ipx dxdt. 

i n=0 3xj 


Now, using Cauchy-Schwartz inequality repeatedly, we obtain 




N-l t"+i 

pX 1 ^ V 


= 


/ ' ^ {r]ik^+i,u])-q{k^_i,u])j(pt dxdt 



j 

J Xj 


j ' 




0 ^ x-j 


^3 + k ^3-h 


\ipt\ dxdt 


^4_L 1 - ^'7— ^ 

J' O. J ' 


(Ax)^ 


l^tl" dx ) dt 


<C{Ax)"^ J j 5 ]|A:^+i j lj2j'^^\V’tfdx\ dt 


<C{Ax)i \k\^y U /2 

A similar argument can be used to show almost verbatim 

i^AliP’d),^) < C(Ax) 5 ll‘/5|l//i(nT) ■ 













14 


R. DUTTA, U. KOLEY, AND D. RAY 


Therefore, summing up, we have 

I I ^ C{Ax)^ l^BV II V^ll iji (Et) ' 


Thus, 


{£Ax(v,<l)}^^yo compact in iJjJ(nT). 


□ 

Next, we need to show the compactness of 9)} ax>o' First note that we can express 

(5.2) as 

A-l 


{£Axiv, q), <p) =^ Dlvikj-^,u]) (p’^+^{x) dx 

„ N-1 t-+i 

/ d{kj-i,u°)ip°{x)dx-Ax'^'^ D_q{kj_i,u]) dt, 

j "'■fj j n=0 

where ipj_i{t) := (p{xj_i,t), and ip‘^^^{x) := (^(a;,This further implies that 

N—1 px - 1 

-{^AxiTi.q),v) = EE r ^ (D\T]{k^_i,u^) + D-q{kj_i,u^)\ ip dxdt 

,u”) {p^_i(t) - p{x,t)^ dxdt 

iV-i + ^ ^x 1 

+ EE ^ D*^_r]{kj_i,Uj) [p’^~^^{x) - p{x,t)) dxdt 

+ E / Vikj-^,u°)p°{x) dx 
■= q)^^) + q)^^) 


N-l t-+i 

+EE 

j n=0 

N-l ^t"+i ^x._^i 

:e 

j n=0 


N-l ^t"+i .x.^i 


/a; 1 

^ 2 


(^D\r]{k^_i,u'^) + D-q{kj_i,uJ)j p dxdt 


where 

(5-5) (f (’7, 9), ‘zs) = XI E 

i "=o 

while {S]^{r], q), p) corresponds to the three remaining summation term. The compactness of 
and {^A^id, q )]will ensure the compactness of '?)}Aa;>o 

Lemma 5.3. 

(5.6) \^Ali'qiq)\ is compact in H^I{Ut). 

< J ArE>0 

Proof. Firstly, we split the SA^iVjq) as 

{^Alid,q),p) ■■= {£A^'^ir],q),y^) + {£Ax‘^ip,q),p) + ^A^'^{d,q),p) 


{^Ax^ip^q)^p) = J2J2 


N-l t"+i 

- ' / 2+f 


!)_ 


j n=0 ■' •- 

N-l ^t"+i fXj+l 

{£Ax^iv,q),p) ='^Y1 / / ^ Dlp{kj_i,u]) (p'^+^ix) - p{x,t)) dxdt 

{^Ax^ip^q)^p) = E / v{kj_i,u°)p°ix) dx. 

x I i 


q{kj_i,Uj) (pj_i{t) - p{x,t)j dxdt 


where 
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We first estimate the term as follows: 


j,n 

A-1 4-+1 

= EE/ 

3 n= 0 “'‘" 


Ax 


— ip{x,t)^ dxdt 


+ 1 


1 - 5 ’ 3- 


Ax 


— ip{x,t)^ dxdt 


EE 

j n=0 


<?),¥>) 

r^+h qikj-i,u^-i)-q{kj_i,u]_i) 


Ax 


— ip{x,t)^ dxdt. 


Since q' and Lp are both bounded, and k is of bounded variation, we see that 

I 9 ))‘Z^) I ^ C'(T) ||<p||iDo(nT) l^lsv- 

Hence we conclude that 

(5.7) IlAx{r],q) G A^ioc(nT)- 

Now the term lAxiv^ q) can be estimated, using Cauchy-Schwartz inequality repeatedly, as follows 


A^-l rt-+^ 

I O 


\{lAx{r],q),p>)\ <CJ2Y1 

3 n= 0 -"*'‘ •' 

_ N-1 /-x I 

<(A.)icEE/ 

A-l ^t"+i 

<{Ax)icY,Y. 

j n=0'''^‘‘ 

( A-1 

At EE 

n=0 J • 


^-^-1 


Ax 


^-^-1 


Ax 




\ipx{d) \ dOdxdt 

5-i 

\ 1/2 

|(/Ja;( 0 )|^ I dxdt 
1/2 

I 1 ^ da: I dt 


dt\ WipW^i 


_f/i(nT) ■ 


Therefore, in view of the a priori bound (4.3b), we reach to the conclusion that 
(5.8) { 7 Aa;(??, 9 )}Aa ;>0 IS compact in (Ht). 

To sum up, making use of (5.7) and (5.8), along with the Lemma 2.2, we conclude that 


(5.9) 


{^AxHv^q)]^ is compact in 

J Ax>0 


A similar type of argument, with the use of Cauchy-Schwartz inequality, yields 

1/2 


N 


(4E(d,g),¥>) <C{AtY/^\AtAx^Y.T.(^>l)'\ llT’llz/bnC- 

\ n=0 j 

Making use of the a priori bound (4.3c), we conclude that 


(5.10) 


{^Ax’^id, d)} is compact in 

^ J Ax>0 


Moreover, note that the bound 

ensures that 

(5.11) 


{£Ax^id,q),‘p) < c* Ill’ll l=o(r) , 
{4E(d,g)l GMioc(nT). 

^ J Ax>0 
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Using (5.9),(5.10) and (5.11) along the lines of Lemma 2.2, proves that < > is com- 

pact. □ 

Next we estimate the remaining term We first note that Taylor’s series expansion 

yields 

(5.12) r]uikj_i,u])Dlu] = DXq{kj_i,u]) - ,0”^")(i:>+^^”)^ 


for some 9^'^^ between u" and At this point, we shall make use of the fully-discrete scheme 


n+1 


{D\'q{kj_i,u^) + D-q{kj_i, f dxdt 


(3.2) and (5.12) to decompose the term £^„,{ri,q) as follows: 

/'“i+i 

EE/ /^^' 

where 

N-l , 


EE 


j n=0 

N-1 


+ D_q{kj_i,u])^ ip dxdt, 

(^Aa; ’^(??,'?),V5) = P r]u{kj_i,u^)Ax {D+D-v^) ip dxdt, 

j n=o "'L 

jV—1 

(^Ax = 7 EE/ / r]u{kj_i,p)fi{Ax) {D\D+D-u^) ip dxdt, 

j n=o 
AT—1 +TI+1 

-'J_/ /** r /\f 1 

(^Ax I -YVuu{kj_i,9p^){Dlujf ipdxdt. 


j n=0 ' 


Our aim is to estimate each of the above terms suitably. In order to do so, we introduce the 
notation := ip(x,t)dx. 

Lemma 5.4. 


(5.13) 


1)1 is compact in 

^ J At>0 


Proof. We rewrite {q, q) as 

{^Ax^i'n,q),‘p) = EE/ 

J n=0 
N-1 

= E/ E 

n=0"'‘" j 


N-1 

= E/ E 

n=0 i 


-q^{k^_.,p)D.{k^^rJpi) + D.q{k^_.,u]) 


dt 


qik^_i,u]) - q{kj_3,u]_^) 

^ Sc 

-riu{k^-l,p)k^_.Djp. 

qik^_i,u]) - q{k^_i,u]_^) 


dt 


Ax 


dt 


N-l ^*’*+1 

e/ e 


n—0 ' 
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Ax 

Since r/' and q' are bounded, making use the total variation bound for k, we conclude that 

<C'|A:|sy IIV5|lL=o(n,j,) • 


<i)j dt 


Hence this implies 
(5.14) 




Next we move on to estimate (rj, q). At this point we shall make use of specific structure 

of q and q given by (2.1). 

Estimate for {q,q) := {qi,qi): First, we recall that 

qi(u) =u-c, qi{k,u) = f{k,u) - f{k,c) = k{x){f{u) - f{c)). 

Thus, we find that 

d ,, (7i(fc^_i,Mp - <7i(fcj_i,u^_i) 




Ax 


fj+L - f^-L / < - / u" 1 

^2 Ax -^2 Ax 

(/(.")-Aj) - (/(«"_,)-/Uj) 


Ax 


= k iD_ 

J 2 


(/(“?)-Ai 


We insert this in the expression of (^Aa;’ ’ (^i9):¥’)' Then, using summation by parts we obtain 

^'-1 ^t"+^ 

dt 


(4Td>n.<n),rf = E / («“")- Ai) 

n=0 j 

W-l 

E I E (/("?) - Ad d.-Ad 


dt 


n—0 
N-1 .fn + l 


E/., E (/("?)-Aj) 


n =0 ' 


kj+i D+kj_i 


dt 


■■= {iAx{qi,qi),A) + {iiAxdni,qi),A)- 
To proceed further, we first estimate the term I/dAxi'Hii Q.i) follows 

N-l 

■ E / E (a“?) - Ad ‘..u ° 


\{lAxivi,Qi),‘P)\ = 


+ dt 


n —0 


< 


'Iloo E / EI« 7 -“”+iI 


n=0 ' 


N-l 

'iloo E / EK”-<-il 


n=0 ' 


Next, we estimate as follows 


- $,-1 = 


r^+i r^-h r^+h / n r^-h ( \ 

/ yjrfx - / (pdx= dx - / 1 dx 

J X . I j X . ^ j X . I ^ ' j X . 3 ^ ' 


X .,1 nx 


(fix dO dx — 


ifx dO dx. 
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Therefore, 


rx.,1 px 


1$,< 


\(px\ dOdx 


nX . 1 nX 


<2 Ax 


\lPx\ dx <2 (Ax) 2 


\^Px\ dO dx 


\(px\ dx 


1/2 


Thus, we conclude that 


(5.15) 


< 2(Ax)1/2 [ [ \^x 

V^^-h 


1/2 


dx 


Therefore, applying Cauchy-Schwartz inequality and using the a priori estimate (4.3b), we obtain 


N 


1/2 


|(/Aa;(T?l,gi),/5)| < 2 ||fc||^ ll/'ll^ (Ax) 1/2 [ ^ AtAx||£i_u" 
< C(Ax)^/2 M\m{nT) ■ 


\H^{U7 


\n—0 


Hence 

(5.16) {1^Ax(??i,9i)}ax>o is compact in 7J"^^(nT). 

On the other hand, making use of the total variation bound for k, we have 

\{IlAx{r]i,Qi),^)\ < C ||v5||ioo(n,r) \^\bv ■ 

Therefore 

(5.17) {77A.(?7i,gi)}A.>0 e Mioc(nT). 

Using (5.16), (5.17) in tandem with Lemma 2.2, we get the compactness of |£a)e’^’^(» 7 , 9 )| 

in H{;l(nT), for (ry,g) := ( 771 , 91 ). 

Estimate of for ( 77 , 9 ) := ( 772 , 92 ): We recall that 


Aai>0 


hence, 


r] 2 {k,u) = k{f{u) - f{c)), 92 (fc,u)=^ 


^m{k,u) = kf{u), ^ 92 (A:,u) = k'^f'{uf. 


This implies that 
d_ 

du 


^ [92(fc,_i,0-92(fc,_i, 77 )^ 1 ) 


1 




Again, by Taylor’s Theorem 

92(fcj-i,u”) -92 (A:j-1,m"_i) 




1 92 

1 52 


= fc?_i/'(<)^ (< - <-i) + ^T2'Z2(fc, 4 ,07 0 « - nUf 


= kun^i) (/«) - /«-i)) + hi^92(fc,_i,r i)«- uu)\ 
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where 0”_i lies between and m"_i. Thus 

J 2 J J 






Therefore 


A-l 


(4’^’'’'(ry2,g2),^) = E / E (fc|-i/'K)^-(/K) - /JVi)) 

n=0 1 


dt 


<Cax('/2,92).<^) 

^-1 + ^ 1 52 


H / -'‘Ur 


n—0 ' 


(Ci^(i;2.92),v> 

Making use of the a priori estimate(4.3b), we conclude that 


A-l „t"+i 


1 92 


(QL(^2,g2),<p) := ^ I E^^«2(V^^^i)K-<-i)' ll<Plloo- 


n=0 


Hence 

(5.18) {QL(^2,g2)}^,>o G Mioc(nT). 

On the other hand, using summation by parts, we can estimate the other term as follows 
v-i ^ 

dt 


(QL(^2,92),^) := E / E (fc|-i/'«)^-(/«) - f^+o) 

n=0 -I*" j 

N-1 t"+i 

= -E / (/K)-/yj)* 

n=0 J 

V-1 t-+i 

= - E / E(/«) - /;+i) iD+k]_,)f{u]^,) $,+1 dt 


n—0 ' 


V-l „t"+i 


{XA^{r]2,q2),v) 


-'Ll E(/W)-4i)bGi (D+/'K))u+.<ii 

n=0"'‘’' .1 


(Vax(»72,92),¥2> 


A-1 


- If / If(/K)-4i)tb. /'K))c-(u)<ii- 


n=0 

"-V-" 

(2ax(i;2,92),V> 

Using essentially the same type of arguments as before, we can deal with the above terms and 
prove the following estimates: 

|('Ta2;(?72,<?2),<^)| < C ||(,9||^ |fc| BV Halloo ’ 

|(3^Ax(r?2,d2),‘/2)l < C \\(p\\^ ||fc||^ , 

\{ZAx{m,q2),^)\ < cM hi ||fc|loo. 
with the help of (5.15), and the a priori estimates (4.3b), and (4.3c). This implies that 
(5.19) {TAa:('^2) 92)}Aa;>0 ’ Ax{V2i ^ 2 )}Ax>0 ^ .^loc(n 7 ’). 
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and 

(5.20) {2Ax{m,q2)}Ax>0 is compact in 

We use Lemma 2.2 along with (5.18), (5.19) and (5.20) to conclude the compactness of 9)| 

in -f^ioc(nT), for (ry,g) := ( 772 , 92 ) as well. 

Finally, the compactness of < 9 ) 1 along with (5.14) ensures the compactness of 

I J Aai>0 

Next, we tackle the term £Al.''^{r], 9 ). 

Lemma 5.5. 


Aa;>0 


(5.21) 


{^Ax^iv,q)]^^ is compact in 


Proof. We first recall the expression for {£ax’ 

N-l „t-+i 


Pu{kj_i,n!^)lAx{D+D-u^) ipdxdt. 

_n I'i 


j n=0 

Using summation-by-parts formula, we can write 

Af-l „t"+i 


{£Ax^iv,q),p) = -^ 


N-l .("+1 


E /„ d-<^ 3 dt. 


n—O ' 




Now we write £Ax’'^'^('d: 9 ) as 


A-l .("+1 


(^Ax^’^(^> 9)><^) =-/3 E / Aa;^U_7 


n=0 


^7?«(fcj--l,Up -77„(fcj_l,M"_l) 

Ax 


dt 


N-l t"+i 

-/^E/ E^^^- 


n—O ' 


(Ai\:c(ri,q),v) 

riu{kj_i, u^_ 1 3 , 7 i”_ 1) 

Ax 


dt . 


We estimate the terms AAxiv,q) and SAa;(77,9), using a priori estimate (4.3b), as 

Af-l .("+1 


\{BAx{v,q),p)\< ^ [ E^^l^- 

n=0 j 


'N-l „t"+i 


— fc,_3 

3 2 J 2 


1/2 


dt 


1/2 


< fd ii^^iioo E / E(^^)' E 

\ n=0 1 / \ 1 


< C/3 ||97||^ \k\sy. 
and, similarly 

|(.AAa,(77,9),</5)l < /? 


N-l .("+1 


AaiE 


|U_u"f dt</3C ||(/.||, 


n—0 ' 
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Hence 

(5.22) 

'?)}Aa:>0 ’ '?)}a£c> 0 ^ •^loc(nT) 

’ 1 , 1 , 2,2 
'Ax 

(4.3b) and the estimate (5.15), we obtain 

N-l 


J Ax>0 

Next consider {rj^q). Using Cauchy-schwartz inequality along with the a priori estimate 

1 the estimate (5.15), we obtain 

N-l ^ 

-n ^ 


n—0 


N-l t-+i 

e Ax\D-u'^\ |Il_<l>j| dt 


n—0 


N-l „t"+i 


1/2 


< pc (Acc)^/^ ( Ax ^ 

\ n—0 ^ 

<pC{Ax)^^‘^ ||(/?||^i(n,j,). 


|D_w"ir dt 


II<P|| 


m(nT) 


Therefore 

(5.23) 


Aa;>0 


{^Ax^’^iV:Q)] ^ is compact in 

^ J Ax>0 

Thus, (5.22), and (5.23) ensure the compactness of |^a)e’^(i?, (z)| 

Next, we focus on £^’^’^( 77 , 9 ). 

Lemma 5.6. 

(5.24) 


□ 


compact in 

Proof. Recall that 

Af-l p 

{£ax ^) = 7EE / r]u{kj_i,u'^)p{Ax) (Zl(i_Il+Il_up ipdxdt. 

j n=o P 

Making use of the summation-by-parts formula and discrete Libnitz rule, we rewrite 

{£ax^{v,(i),p) = -7 E / D_ l-^vikj_i,u]) j dt 

n=0 •'*" ’■ \ “ / 


(Cax(i;,9),v> 


V-l t"+- ^ 

n=0 I “ 


{X>Ax(T?,g),¥’) 

Again, making use of the a priori estimate (4.3d) reveals that the term X>ax(^, can be estimated 


as 


/ N-l \ 

|(2?Ax(i?,g),/3)| < 7 C (Ax)5 I fi{Ax)AxAt I llv^ll^/qn^r) - 7C'(Aa;)5 ||¥>||^i(n,r) 

V n=0 / 

SO that 

(5.25) {d^Ax(i?,g)}Aa :>0 IS compact in i7j;;](n'r). 

On the other hand, to estimate C\x(jIj q) term, we first split that term as 

{Cax{v,q),p) 
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= -7^/ ^x{l^x)D\D_i 


Ax 


dt 


{G&n:{v,g),‘p) 


- , Y. c m(ax) o^D^uf-AizklAA^AhziAA. t.._, , 


(^Ax(i;,9),<p> 

A similar type of argument, as used before, can be used to deal with the terms and J^Aa;) with 
the help of the a priori estimate (4.3d), returns 

\{Q/^x{'n,q),^)\ 

( "" 

<27 ||<p||ioo(nj,) ( X! At Ax 

and 


n—0 


1/2 


N 


1/2 


'^AtAx\\D_u^f \ < 7 C ||<p||ioo(n,j,), 


\n—0 


N-1 


1/2 


|(-^Ax(??,9),</5)| < 7BV(fc) [tJ-i^x)AxAtJ2\\DXD_u'^\\ < 7 C 


n—0 


This implies that 

(5.26) {^Aa:(7; 9 )}ax>0 > 9)}ax>0 ^ Iq;. (ny) . 

Therefore, in view of the Lemma 2.2 along with (5.25), and (5.26), we reach at the conclusion 


(5.27) 


is compact in 


□ 


Making use of Lemma 5.4,5.5 and 5.6, we can finally prove the following. 

Lemma 5.7. 


(5.28) 


{■S’ax(^:9)} 


Aai>0 


is compact in 


Proof. We have already shown that |£^'^’^(? 7 , g) I A^Ax^{'qyQ)\ and |£’^’j;’^(? 7 , g)| 

L J Aic>0 ^ J Ax>0 L J 


are compact H ^(LIt). Thus, we only need to find an estimate for £as’^(7j '?)■ using the a priori 
bound (4.3b), we have 


Aai>0 


N 


{^Ax^iv,q),p) <C'll<7’lloo ^cIIt’IIoo^ 


\n—0 


so that 
(5.29) 


{4’^’"(7,g)}, GMoc(nT). 

^ J Ax>0 


An application of Lemma 2.2 allows us to conclude (5.28), thus proving the lemma. 


□ 


Combining Lemma 5.3 and 5.7 ensures that 
(5.30) {^Ax{v,q)}^^yo is compact in iJ-i(nT). 

Thus, Lemma 5.2 and (5.30) justify the assumptions we had made in the proof of Lemma 5.1. 

Now we are in a position to state the “convergence theorem” which guarantees the convergence 
of approximate solutions {?iAa:}Ax> 0 ’ generated by the scheme (3.2), to a weak solution of (1.4). 
The following theorem can also be viewed as a modified version of the classical Lax-Wendroff 
theorem (for more details, consult the monograph by LeVeque [21]). 
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Theorem 5.1. Let u^x be a sequence of approximations generated via the seheme (3.2) with 
fi{Ax) = 0[Ax^). Then there exists a funetion u S L°°([0, T]; L[qj,(R)) and a subsequence of 
{Ax} (not relabeled) such that uax u as Ax } 0. Moreover, the function u is a weak solution 
to (1.4). 


Proof. The strong convergence of uax to a function u immediately follows from the Theorem 2.1 
and the Lemma 5.1. It remains to show that m is a weak solution. Let ip G C^(R x [0,T)) be 
any test function and denote ip" = ip{xj,t"). Multiplying the scheme (3.2) by AxAtip", and 
subsequently suming over all j and n yields 


AxAt 




1+5 


= P AxAt ^ ^ Ax Ip" Dj^D_u" + 7 AxAt ^ ^ /r(Ax) ip" D*^Dj^D_u". 

j n j n 

By standard arguments, it is clear that 

AxAt ^ ^ ip^ DXu^ = -AxAt Y: E E € 




- 


/ / uiptdxdt— / '0(x, 0)tto(2;) dx, as Ax } 0. 

JM. Jo Jr 


Next, using the a priori bound (4.3b), we conclude 

P AxAt E E D+D_u" = —P AxAt E E 


3 " 


3 n 


1/2 


1/2 


< P AxAt y/ A,Ty/(D_,.“)d \AxY,(v-t 


n\2 
3 


<P Ax^AtYY^^- 


n J 


1/2 / ^ 1/2 
■"'2 1 i Ax^AtYY^^-^3^^ 

n j 


< /3C'(Ax)1/2 ||i/;||^2((0,T);ffi(R)) ^ 0. asAx}0. 

Repeated use of Cauchy-Schwartz inequality along with the help of a priori bound (4.3d) returns 
7 AxAt E Ed(^^) D(^_D+D-u" 

3 " 

1/2 


< iC 


/l(Ax)^/2 

Ax1/2 


E] E/ M(^^)^2;2At \d(^D-i 




m(Ax)^/2 

^^ 1/2 ll'*/'llL 2 (( 0 ,T);ffi(R)) 0; as Ax I 0. 

Finally, a simple use of summation-by-parts formula implies 

AiAt^^VJD- (7+p;h) = 


3 


3 ^ 


= - AtE ^ D+'ip] k{x) f{uj) dx - AtE ^ D+ip)) (kj^i - k{x)^ f"_^_i dx 


3 ,n 1 


J," 1 


- At 


r, ' 1^1 


3 kl 1 


1/2 
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where, again, a standard argument reveals that 

oT 

■'Ax 




/ / k{x) f{u)'t{jxdxdt, a,s Ax 10. 

Jr Jo 

Y. I m * < f IM.T) - *+^ /■ - *<">! 

„ Jii ^ Jr ^ Jr 


Now observe that 


Aa:/2 


dx 


and consequently 


/ k{x)-k 


d+2 


dx !->■ 0 as Ax I 0. 


Keeping this in mind, we find that 


f-i — 
^Ax — 




j,n 


< M 


Li([0.T];L° 


H / k{x) - fc^+l 


dx !->■ 0 as Ax I 0. 


Finally, making use of the a priori bound (4.3b), the last term can be estimated as follows: 


fi — 
^Ax — 


D+'<P]'k{x) -/(mj)) dx 


j n 


1/2 


< M(Ax)^/2 I 


lL2([o.T];ffi(B)) 


Ax2At^^(d? _u")^ j H> Oas Ax I 0. 


3 " 


To sum up, we have proved that u is a weak solution of (1-4), i.e.. 


f f {tptu + iljxk{x)f{u)) dx dt + f i/>(x, O)uo(x) dx = 0, for all i/> S X [0, T)). 

Jr Jo Jr 


□ 


6 . Entropy Solution 

As we have already mentioned in Remark 3.2, we use the Engquist-Osher scheme to make the 
analysis more concrete, but our methods can easily be adapted to general monotone schemes. 
Drawing preliminary motivation from Karlsen et al. [14] and Towers [28], we first show that the 
limit solution satisfies the Kruzkov entropy inequalities locally, away from the jumps in k. Then 
we proceed to show that the limit solution satisfies Kruzkov type entropy inequalities when the 
test function has support which intersects one or more jumps in k, and finally we show that this 
implies Kruzkov entropy solution. 

To proceed further, we need some additional regularity assumptions on the discontinuous coef¬ 
ficient fc(x). We assume that k{x) is piecewise Lipschitz continuous in K with finitely many jumps 
(in k and k'), located at ■^ 1 )^ 2 ) • • • iCm- More specifically, we assume that there are finitely many 
Lipschitz continuous curves a;i,a; 2 , • • • ,u!m such that G oji, the union of which we denote by 

M 

Wj — D. 

i=l 

For the sake of simplicity, we also assume that none of the curves intersect. The curves wi, 022 , • • • , ujm 
partition K \ D into a finite union of open sets: 

K \ n = 'R.q U TZi U ■ • • 'R-Mj 

with the curve separating the sets 'R-m-i and TZm- We will assume that 

k G Lip(7?,„), m = 1, 2, • • • ,M. 
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With this assumption, k has well dehned limits from the right and left along each of the curves 
(jjm, 1 <m < M, and we denote these limits by := k{^^) respectively. 

To this end, we first establish the following entropy inequality for smooth entropy function rj, 
to avoid complications arising from the discontinuity 77 ' (u) = sign {u — c) in Kruzkov entropy. 


Lemma 6.1. Let{r],Q) be a convex entropy pair with rj being a -function and rj'(u) =Q'{u)f'(u). 
For the test function ip > 0 with compact support int > 0, x G K\r2, and every c € K, the following 
entropy inequality holds 

/ / (r]{u)ipt-\-kQ{u)ipx)dxdt-\- r]{uo)ip{x,0) dx 
JwJo ^ ' JR 

~ f f k'(x) (j]'(u) f{u) — Q{u)^ip dx dt > 0. 
Proof. We begin by rewriting the scheme (3.2) as 


= wf +/SAxAtU+U-u" + 7 At/i(Aa;)ilj 


with 


7 = 7-AA+(t,_j/jLj) 


where A+ denote the undivided forward difference operator, and /", 1 is the EO flux corresponding 

(?+ 2 

to the flux function f{u). Let FIj^i be a entropy flux function, consistent with Q{u). Then, in 
light of Kenneth et al. [14, Lemma 4.1], it is evident that for any c S K 

K - c| < |u" - c| - A i 1 ) - Asign {w^ - c) /(c) A+fc^_i. 

To obtain a similar type inequality for a smooth entropy-entropy flux pair, we follow a classical 
approximation argument. For a rigorous proof, we refer to the paper by Karlsen et al. [15, Lemma 
5.7]. In what follows, we have the following inequality: 

Viwj) < !?(«”) - A (fcj+i ilj+i - kj_iHj_i^ + A (Q(w'*) - g'{w^)f{w^)) A+fc^_i 

After rearranging terms, a discrete entropy inequality for the scheme results 

< 7 + - fK)) - AAh- + Af“A+i,_.. 

where = Q{wf) — il'{wf)f{wf) and 7 " = Next, multiplying the above inequality by 

AxAtipf with ipf = ip{xj, tn), where ip smooth, non-negative test function ip with compact support 
in (K \ fl) X [0, T) and using summation by parts yields 

AxAt EE Ip" 7?+r?” - AxAt 

j n j n 

- ^^ ip- + (lyj+i - 7(0) /At) < 0. 

j n 

As before, a standard argument reveals that 

AxAt EE iP// DXg{uJ) = -AxAt Y. Y iK) - i: 7“ ’)(«“) 

j n j ri j 

!—>■—/ / r]{u)'iljt dx dt — / 0 )? 7 (uo(x)) dx, as Ax st 0. 

Jr Jq Jr 

Next, we rewrite 

= Aa:At^^Oi (Oj -QK)) 
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+ Ax At 


where it is straightforward to conclude that 

cT 




/ / k{x)Q{u)'4!xdxdt,aaAxiQ, 
Jts. Jo 


and the a priori estimate (4.3b) gives that 

I^LI <C(Ax)i/2 (Ax^AtEl^- 

\ n,j 


1/2 


Il2([o,t];H1(r)) 0 as Ax I 0. 


To estimate the term involving D+kj_i, we split the term as follows 
AxAt = AxAtEE (QK) - ^'K)/K)) 

n j 

= AxAt E E (^?K) - ^'K)/K)) 


" 3 


AxAt EE (Q«) - Q«)) 


+ AxAtEE (^'K)/K) - v'{u])f{u])) . 


Again, a straightforward argument shows that 


^Ax 


[ [ {Q{u) — r]'{u)f{u))k'{x)'ijjdxdt, as Ax ), 0. 
Js. Jo 


For the other terms, we intend to show that 

^Ax^^Ax 0 as Ax i 0. 

To see this, first note that 


(6.1) 


IQ)^”) - Q{vj;)\ < C U” - u”! < CA 


< CA 




k A k A _1. 

2 J 2 


■^i+1 “'i I ri-i 

Since k{x) is Lipschitz continuous within the support of '0, it is clear tha 




A:,-i 1 — fc,_ 1 

JV 2 J ; 


D+k 1 

^ J 1 


< 


^Ax||fc'||^ \k\j^y hA 0 as Ax 4.0, 


and 

AxAt EE^" 


dd+k,_i 


-‘i+i “i 


n J 


1/2 


< (Ax)i/^ Ax^AtEl^-<n ll^'lloo IIV'IIl2([0,t];l2(h)) ^ 0 as AxiO. 


This proves that i—)■ 0 as Ax i-A 0. A similar calculations show that i-A 0 as Ax i-A 0. 
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Finally, the last term can be estimated via Taylor series as follows 

AxAt EE ^7 


n j 


= P AxAt ^ ^ 'tl;^7]'{u^)AxD^D-u'^ + 7 AxAt ^ ^ {uJ)fi{Ax)D:^D^D^u'^ 

j n j n 

- axY,xv"{0,) 7" A+ (7-i/7i) (7^' - O + ^^^E - O' 

ii’T' 

:= Qax + 2ax + Qax + Qax- 

For the first term Q\^, using the discrete chain rule, we proceed as follows 

QX, := /? AxAt E 7”^'(0 ^^D+D_u] = /? AxAt ^7 (^'(O ^-^7) 


j n 


j n 


- /3 AxAt ^ ^ Ax 7" D+u] 7(7) ^+7- 


j n 


where using the non-negativity of the test function -i/j, we conclude that 


4 . > 0 . 


Moreover, using the a priori bound (4.3b), we find 


\£ax\< PC{Axf/^ ^Ax^AtY^'^{D_u]f"^ |772([0,t];J71(K)) ^ 

For the next term we proceed as follows: 

Qax ■= 1 ^xAtY^^-ipJr)' {u]) ^i{Ax)D7 D+D_u] = 7 AxAt /r( Ax) 7"' D+ [r]' {u^) D+D _u^) 


3 " 


el. 


— 7 AxAt y: e 7” DtD+u] 7(7) ^+7 • 


fix 


Making use of the a priori bound (4.3d), we conclude 


\£lx\<ic 


1/2 


fi{Ax) 

Ax^/^ 


1/2 


^(Ax) Ax^ At 


,n\2 


llL 2 ([ 0 ,T];ffi(B)) H> Oas Ax 4, 0, 




and 


I^LI <7C 


/x(Ax)^/^ 


1/2 


1/2 


^Ji{Ax)Ax^AtYY^^^^-^7f Ax^AtYY^^-^7f 


3 " 


Ax 

i-> 0 as Ax 4, 0, 

Next, we turn our focus on the term Q\^. In fact, we write 

ei, := Ax 5 : 5 : A,"(« 2 ) 7” A+ (u_./-_.) (»"+' - 0 

j n 


i n 
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= /3Ax^^A 77"(02)V-”A+ (Vi/T-i) ^^^tD+D-u] 

3 ^ 

'-s/-" 

‘^Ax 

+ 7Ax^5]Ar?"(02)V-,"A+ (Vl/T-i) AiMAi) 

j n 

" -V-^ 

fix 

Before we proceed, we recall that in accordance to Remark 4.2 A = (^{Ax)/{Ax)'^, or in other 
words At = (^{Ax)/Ax. We also need to use the identity AxD^D_u^ = — D-U^- Thus, 

using the discrete chain rule for the term A+ and apply Cauchy-Schwartz inequality 

repeatedly, we obtain 

/ \ 




Ax" 


+ 




J ri 


i-A Oas Ax I 0. 


Similarly, 


|£-« 


n(Ax) 


1/2 


l^lsY Ax^/i(Ax)At^^(D^|.D_i 


.n\2 


. AxV(Ax)A t ^ 'Y^{D\D_u^f 


i-A 0 as Ax I 0. 


Finally, we are left with the term Q%^. We see that 

eL ;= -O” 


J n 


< p‘^Ax EE^” 77"(6Ii) {AxAtD+D_i 


+ 7 ^ Ax EE^^"(^o (^fi{Ax)AtD^ D^D_i 

3 n 

' -V- 

c-10 

^Ax 

We start with the first term Making use of the a priori estimate (4.3b), we conclude 

,m(Ax) 


I^LI <c^ 


Ax^ 


Ax^At^^(L>_upM ||i/;||^ hAOasAxiO. 


Similarly, making use of the a priori estimate (4.3d), we argue that 
I yi 

C A 


?A.| |^AxV(Ax)At^^(Z?^Z?_u^")2j lliAII^^OasAxiO. 


□ 


Now in view of the above result, we are ready to prove the following lemma: 
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Lemma 6.2. Let u{x,t) be a weak solution constructed as the limit of the approximation uax 
generated by the scheme (3.2) with iJ,{Ax) = o{Ax'^). For the test function if > 0 with compact 
support in t > 0, a; G K \ and every c G M, the following entropy inequality holds 

T 

{\u- c\'ipt + sign {u - c) k{x) {f{u) - /(c)) ifx) dx dt 



+ / \uo 

Jm 


c\ ip{x, 0) dx 



sign (u 


c) k'{x) /(c) ifdxdt> 0 . 


Proof. A simple manifestation of the above Lemma 6.1 for the specific entropy r]{u) = \u — c\ and 
consequently the entropy flux function Q{u) = sign (u — c) (/(u) — /(c)) essentially completes the 
proof. □ 


Lemma 6.3. Let u{x,t) be a weak solution constructed as the limit of the approximation uax 
generated by the scheme (3.2) with p,{Ax) = o{Ax^). Let 0 < if £ I?(]R x [0,T]). Then the 
following entropy inequality is satisfied for all c gM. 

/ / {\u-c\ ift + sign {u - c) k{x) {f{u) - f{c)) ifx) dxdt + / \uo - c\if{x,0) dx 
Jr Jo Jr 

+ l/(c)|/ / \kfx)\ifdxdt+Y^ / |/(c)(fc+ - fc;;)|'0(Cm,t)dt > 0. 

JR\n Jo _ 1 ^0 


m—1 


Proof. A straightforward adaptation of [14, Lemma 4.2] along with the help of Lemma 6.1 con¬ 
cludes the proof. □ 


To proceed further, we combine above two lemma’s. In what follows, we have the following 
important theorem: 

Theorem 6.1. Let u{x,t) be a weak solution constructed as the limit of the approximation uax 
generated by the scheme (3.2) with yb{Ax) = o{Ax^). Let 0 < if G I?(K x [0,T]). Then the 
following entropy inequality is satisfied for all c gM. 

/ / {\u - c\ ift + sign {u - c) k{x) {f{u) - f{c)) ifx) dxdt + / \uo - c\if{x,0) dx 

Jr Jo Jr 

j- pT M T 

+ / sign{u-c)k\x) f{c) if dxdt+'^ / |/(c)(fc+- fc;;)|t) df > 0 . 

Proof. A verbatim copy of the proof of [14, Lemma 4.4] ensures the proof. 

□ 


6.1. Uniqueness of Entropy Solutions. Following [14], we mention that entropy solutions 
are unique under a crossing condition. It is well known that one has to impose the crossing 
condition only because the entropy inequality ( 1 . 6 ) alone is not sufficient to guarantee uniqueness 
when the crossing condition is violated. However, we mention that in the multiplicative case 
f{k, u) = k{x) f{u) there is no flux crossing, hence we don’t assume any such crossing condition. 

One more technical issue is the existence of traces along the discontinuity curves Wm, for 
m = 1,2, •• • , M. Without going into the details of it, we simply mention one instance where 
we automatically have the existence of strong traces. We remark that, due to the genuinely non¬ 
linearity assumption, existence of traces is guaranteed when k{x) is constant on each region 
for m = 1, 2, • • • , M. This is a consequence of a general result by Vasseuer [30]. For more general 
case, we encourage the readers to consult [14] for general assumptions ensuring the existence of 
traces. 

To this end, we collect all the above mentioned result in the following theorem. 

Theorem 6.2. Assume that the assumptions A.l, A.2, A.3, A.4, and B.l hold. Moreover, 
suppose that the addition regularity assumptions on the discontinuous coefficient k holds, and the 
traces along the discontinuity curves exist. Then a unique entropy solution to the Cauchy problem 
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(1.4) exists and the entire computed sequence of approximations {uax}ax> 0 ’ generated by the 
scheme (3.2) with ^(Ax) = C>(Ax^), converges to the unique entropy solution of (1.4). 

7. Diffusive-Dispersive Approximations 

As we already pointed out in the introduction, all the aforementioned techniques can be utilized 
to analyze the equation given by the diffusive dispersive approximations of scalar conservation laws 
with a discontinuous flux of the type 


(7.1) 


uf + f{k{x),u^)x = e/3ul^ + xGKx (0 ,r), 

u^(x,0) = uo(x), X G K, 


when e > 0 tends to zero with 0 < /r(e) i—>■ 0 as e i—?> 0. Here T > 0 is fixed, /3,7 > 0 are fixed 
parameters, u® : K x [0, T) i—)■ K is the unknown scalar map, Uq the initial data, fc : K i—>■ K is 
a spatially varying coefficient, and the flux function / : i—)■ K is a sufficiently smooth scalar 

function. 

We propose the following fully-discrete (in space and time) finite difference scheme approximat¬ 
ing the limiting solutions generated by the equation (7.1) 

-t A>-h” 1 = l3AxD+D_u] + jyi{Ax)D+Dlu], j eZ,ne No, 


(7.2) 

(7.3) 


= ^ / " uo{0)d9, j G Z, 


where /3,7 > 0 are fixed parameters, and /i(Ax) i-G- 0 as Ax i-G- 0. As before, we will either use 
/i(Ax) = 0{Ax^) or /r(Ax) = C>(Ax^) depending on the quest for the convergence of approximate 
solution uax towards a weak solution or the entropy solution, respectively. Note that in this case, 
in contrast to the a priori estimates in section 4, we have the following estimates 

Lemma 7.1. LetuAx be a sequence of approximations generated by the scheme (7.2)-(7.3). More¬ 
over, assume that the initial data uq lies in L^(IR). Then the following estimate holds 


(7.4) 


\dX\\u- 


+ 6 


7 Axp,(Ax) 


+ JAx||D_u"|r < C, 


provided At and Ax satisfies the following CFL condition 
(7.5) 


max l2\\ P 


Ax^ 


_ ^ g 7M(Aa:) 

7 ^(Ax) Ax^ 


: 8A ||fc|l^ ll/'ll ^ < min(l - i5,/3 - (5), 5 G (0, min (1,/3)), 


where A = At/Ax and the constant C > Q is independent of Ax. 

In particular, the estimate (7.4) guarantees following space-time estimates: 


(7.6a) 

VnGN, Ax^(u”)2<C', 

j 

(7.6b) 

Ax^AtY,Y.^D_u^f <C, 

j n 

(7.6c) 

AtAx^/i(Ax) ^ ^(D^uP^ < C. 
j 

(7.6d) 

Ax2At5]^(D^u7)2<C, 


Proof. To start with, we multiply the scheme (7.2) by Axu" and subsequently sum over j G Z. 
Then, using summation-by-parts formula and the identity (3.1), we obtain 


IaA/) 
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= —/3Aa;^^ Aa; — jfi{Ax)Ax'^^D_u'JD_{D_u^). 


1 IT 

u-D_u]=^D_iu]r + —iD_u]r 


Using the identity 
(7.7) 

which is very similar to (3.1), we get 


=0 

Furthermore, it has already been shown in the proof of Lemma 4.1 that 

< c 

where C is independent of Ax. Thus, we have the estimate 

(7.8) \dX \Kf < ^ ll^Xf - /3Ax \\D.u-f - 
Noting that AxD+D^ = D-{D+ — D-), we use the scheme (7.2) to get the relation 


IdXu'^W +c 


(7.9) 

Now 


— WDXu'^f < At 

O M ~r M — 




+ l3^AtAx'^ \\D+D_u 

7 ^At^(Aa;''^ 


+ 


Acc^ 


111,2 




At 






= ^E(»^,+jAj +‘^,-sAi - 7 -sAj 

3 

i “i: (k,,, - k,_,y (4j)» + (k,_. f (/y, - /y,)" 

3 

Since EO flux is Lipschitz continuous with Lipschitz constant ||/^||oo (c.f. (3.5)), we conclude that 


sup sup 

3 ^ 


jpn 

U+5 


< ll/lloo 


, and 


£71 Jen 

U+4 Jj-k 


<ll/'lloo(h”-«^l| + 


^ 1+1 


Thus, we get 
(7.10) 


At 




< 2A ||fc|| ll/f \k\^y + 4At \\kf wrf ) 


Using relations (7.8),(7.9) and (7.10) in tandem with the fact that ||I?_(.)"’|| = ||Z1+(.)”||, we get 
the estimate 


1 


(7.11) 


Di\\u^f+{l-2 


(3'^XAx^ ^'yXfi{Ax)\ ^Axyu{Ax) 

jfj.{Ax) Ax^ 


\D^ m " 


+ (/3 - 8A Ilfcf ll/'f) Aa; ||iA_u"f < C 
Choosing A in accordance to (7.5) ensures that for some 6 G (0,min(/3,1)) 

B'^XAx^ jXfi{Ax) 


1 - 2 - 


-4- 


>5, and/3-8A||fcf ll/'f >(5, 


7 ^(Aa:) Aa;^ 
thus leading to the estimate (7.4). 

In order to prove estimates (7.6a), (7.6b) and (7.6c), we multiply the inequality (7.4) by At 
and subsequently sum over all n = 0,1, • • • , A — 1 to reach 

1 ||^JV||2 ^ g iAxAUliAx) ^ 11^2 y„||2 ^ SAxAtY,\\D-u^^ < ^ ||wof + C 
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This essentially finishes the proof of the a priori bounds (7.6a), (7.6b) and (7.6c) . Using these 
estimates and the relation (7.9), a direct computation shows that (7.6d) holds. This completes 
the proof. □ 


Note that, making use of these a priori estimates, an appropriate “convergence theorem” can 
be formulated, which guarantees the convergence of approximate solutions {'i*Aa:}Ax> 0 ’ generated 
by the scheme (7.2)-(7.3), to a weak solution of (1.4). 

Theorem 7.1. Letu^x be a sequence of approximations generated via the scheme (7.2)-(7.3) with 
fj,{Ax) = 0[Ax^). Then there exists a function u G L°°([0,T];Lj^^.)®)) and a sequence {Aa;^} of 
{Ax} such that u^xj ^ u as Axj } 0. Moreover, the function u is a weak solution to (1.4). 

Proof. The proof of this theorem is very much similar to the proof of the Theorem 5.1, except 
the analysis of the terms involving Dj^.D^uf. However, these term can be treated like the term 
in section 5, but we need to use the a priori bound (7.6c) instead of (4.3d). For 
brevity of exposition, we omit the details of the proof. □ 


A similar result, in view of the analysis in section 6 and above priori estimates, can be obtained 
regarding the convergence of {uAs}Aa;>o’ generated by the scheme (7.2)-(7.3), to the unique 
entropy solution of (1.4). 


Theorem 7.2. Let u{x,t) be a weak solution constructed as the limit of the approximation uax 
generated by the scheme (7.2)-(7.3) with /i(Ax) = C>(Ax^). Let 0 < ip G DIM. x [0,T]). Then the 
following entropy inequality is satisfied for all c gM. 


/ / {\u-c\ipt + sign{u-c)k{x){f{u)-f{c))'tpx) dxdt+ \uo-c\ip{x,0)dx 
Jm. Jo Jr 

p pT M T 

+ / sign{u-c)k'{x)f{c)ipdxdt+'^ \f{c){k+- k~)\ip{fm,t) dt > 0. 

jR\n Jo _ 1 Jo 


Proof. This can be achieved using similar arguments used in the proof of Theorem 6.1. □ 


8. Numerical Experiments 

We present a few numerical results to substantiate the results we have shown in the previous 
sections. We consider the capillarity problem approximated by the scheme (3.2), as well as the 
diffusive-dispersive problem approximated by (7.2). Note that while the former is approximated 
by an implicit type scheme, the latter is an explicit-in-time scheme. 

8.1. Capillarity approximation. We consider the flow of two phases in a heterogeneous porous 
medium, in the limit of vanishing dynamic capillary pressure. The model equation is given by 

(8.1) ut + f{k{x), u)x = eP {g{k{x),u)ux)x + At(e)7 {h{k{x),u)uxt)x, x S K x (0, T), 

where g, h : —>■'R are assumed to be smooth functions such that 

a < g{.,.),h{.,.) 

for some constant a > 0. This model has also been considered and numerically analysed in [4], 
[16] and [29]. Note that, we have theoretically shown the convergence for the special case when 
g = h = 1. However, we mention that this cosmetic changes in the equation has no effect on the 
central idea of the paper and a straightforwardly incremental modification of our analysis can be 
adopted to analyze the equation (8.1). We have decide to work with this equation because of the 
availability of results for such equations which help us to compare our results. 
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As done in [4], we choose the various quantititis in the model as follows: 

z'^il - k{x)z°) 


f{k{x),u) = 


zW 

g{k{x),u) = k{x)gi{u) = k{x)- 

A 

h{k(x),u) = k{x) hi(u) = k{x)- 
P{u) = ~ ^ 

where z'^ = v? and z° = — From the physical point of view, u{x, t) and (1 —u(a;, t)) represent 

the water and oil saturations respectively, while k{x) corresponds to the rock permeability. The 
corresponding modified numerical approximation is chosen to be 


-P'{u) 


■w _ 1 _ 


(5i(m”) +5i(m"i)) 


= l3AxD+ I 1 ^ D_u" 


+ ig.{Ax)Dj^ fc _i 


(^iK) +hi(M”i)) 


DXD.u^ 


As stated earlier, we can replace the EO numerical flux with any other monotone flux, with a 
similar covergence analysis following through. For simplicity, we choose the Lax Friedrichs flux 
for the present model. The spatial domain is [0, 2] with an initial solution profile 


uo(x) = 


0.8 

0.2 


for x < 0.25, 
for x > 0.25 


Furthermore, we choose /3 = 6 , 7 = 36 and CFL 0.3, with the final time of simulation being 
T = 0.6. 


8.1.1. Continuous flux. We first consider the scenario when fc = 1. This corresponds to the flow 
in a homogeneous rock structure. When the dispersion coefficient is chosen as ii{Ax) = Ax^, 
the numerical solution converges to a weak solution of (8.1), as shown in Figure 8.1. The weak 
solution consists of a leading classical shock wave and a trailing non-classical shock wave, with an 
intermediate state in between. However, when the dispersion coefficient is chosen to fj,{Ax) = Ax^, 
the solution converges to the entropy solution, as shown in Figure 8.2. This is in accordance with 
Theorem 6.2. 


8.1.2. Discontinuous flux. We next consider the scenario depicting the flow through a heteroge¬ 
neous medium. We chose the rock permeability as 


fc(x) 


1.1 for X < 0.6 

1.4 for X > 0.6 


which corresponds to two rock types with a sharp interface at x = 0.6. As before we first consider 
the solution by setting ii(Ax) = Ax^, which is shown in Figure 8.3. The numerical approximation 
is a weak solution consisting of a leading classical shock wave and a trailing non-classical shock 
wave separated by an intermediate state, and discontinuity at x = 0.6 corresponding to rock 
structure. Once again,the non-classical shock disappears when the dispersion coefficient is chosen 
as g{Ax) = Ax^, with the solution approximating the entropy solution. 
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Figure 8.1. Two phase flow through a homogeous medium {k{x) = 1) at time 
T = 0.6 with fj,{Ax) = Ax^. Mesh refinement study indicates convergence to a 
non-classical weak solution of the underlying conservation law. 



Figure 8.2. Two phase flow through a homogeous medium {k{x) = 1) at time 
T = 0.6. A non-classical weak solution is obtained for ii{Ax) = Ax'^, while the 
unqiue entropy solution is obtained for ii[Ax) = Ax^. 


8.2. Diffusive-dispersive model. As we have already mentioned, the convergence analysis for 
the diffusive-dispersive equation (7.1) is very similar to the other. We work with the flux function 

f{k{x),u) = k{x){u^ — u) 

and the numerical scheme given by (7.3). The continuous flux version of this problem has been 
studied numerically in [2] as well. The spatial domain is [—0.5,0.5] with the Hnal simulation time 
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Figure 8.3. Two phase flow through a hetergeneous medium at time T = 0.6. 
A non-classical weak solution is obtained for fj,{Ax) = Ax'^, while the unqiue 
entropy solution is obtained for ii[Ax) = Ax^. 


being T = 0.01. The initial profile of the solution for this problem is taken to be 

f 4 for a; < 0, 

1—2 for a: > 0 

with the model parameters chosen as j5 = b and 7 = 20. For this problem, the EO numerical flux 
is used. 


8.2.1. Continuous flux. We first work with a continuous flux by choosing k{x) = 1. For fi{Ax) = 
Aa:^, the numerical solution converges to a weak solution (see Figure 8.4) consisting of trailing 
classical shock wave and a leading non-classical shock wave separated by an intermediate state. 
Figure 8.5 shows that for the choice /i(Aa;) = Aa:^'®, we get an approximation to the unique 
entropy solution, which corresponds to a single shock wave satisfying Oleinik’s entropy condition 
[24]. 


8.2.2. Discontinuous flux. We work with a discontinuous flux characterised by 


k{x) 


1.1 for X < 0.1 

0.9 for X > 0.1 


The numerical approximation shown in Figure 8.6 with ii[Ax) = Ax'^, is a weak solution con¬ 
sisting of a trailing classical shock wave and a leading non-classical shock wave separated by an 
intermediate state. In addition, there is a standing discontinuity at a; = 0.1 corresponding to 
discontinuiuty in the flux. The non-classical shock disappears when the dispersion coefficient is 
chosen as ^{Ax) = Ax^'^, with the solution approximating the entropy solution. 
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